!SUBROUTINE dynamical
!  USE lapack
!    USE global
!  USE constants, ONLY : c_pi
!  USE forceconstants, ONLY : x_harm
!
!  IMPLICIT NONE
!
!  REAL*8, ALLOCATABLE ::f(:,:,:,:,:)     ! Force constant matrix
!  REAL*8, ALLOCATABLE ::freq(:)
!  COMPLEX*16 ,ALLOCATABLE ::  nc(:,:) ! Normal mode coordinates
!  COMPLEX*16 ,ALLOCATABLE ::  df(:,:),eigenvectors(:,:) ! Dynamical (k-dependent) force constant matrix
!  COMPLEX*16 :: phase
!  REAL*8 :: realk,x,y,z
!  INTEGER k
!  INTEGER cell1,cell2,cell3
!  INTEGER kx,ky,kz
!  INTEGER i,j,m
!  CHARACTER::massfile*80,hessfile*80
!  NAMELIST /param/maxcell,maxk,ndim,massfile,hessfile
!
!
!  READ(5,param,END=10)
!  10 CONTINUE
!
!  !---- Read atomic masses
!  OPEN(50,file=TRIM(massfile),form='formatted')
!  REWIND(50)
!  READ(50,*) kx,ky,kz
!  WRITE(*,*) 'kx,ky,kz = ',kx,ky,kz
!  READ(50,*) Natom
!  WRITE(*,*) 'Natom = ',Natom
!
!  Ndof=ndim*Natom
!  ALLOCATE(nc(Ndof,Ndof),df(Ndof,Ndof),mass(Natom),freq(Ndof))
!  ALLOCATE(f(Ndof,Ndof,-maxcell:maxcell,-maxcell:maxcell,-maxcell:maxcell))
!  ALLOCATE(eigenvectors(Ndof,Ndof))
!  ALLOCATE(k_freq(-maxk:maxk,Ndof),k_coef(-maxk:maxk,Ndof,Ndof))
!  ALLOCATE(x_harm(Ndof,Ndof,-maxcell:maxcell))
!
!  DO i = 1, Natom
!     READ(50,*) mass(i)
!     WRITE(*,*) i,' mass / amu = ',mass(i)
!  ENDDO
!  CLOSE(50)
!  f=0.0d0
!
! !Read force constants from hess file
!  OPEN(50,file=TRIM(hessfile),form='formatted')
!  REWIND(50)
!100 READ(50,*,END=200) cell1,cell2,cell3
!  WRITE(*,*) cell1,cell2,cell3
!  DO i = 1, Ndof
!     DO j = 1, Ndof
!        READ(50,*) f(j,i,cell1,cell2,cell3)
!        x_harm(j,i,cell2)=f(j,i,cell1,cell2,cell3)
!        f(j,i,cell1,cell2,cell3) = f(j,i,cell1,cell2,cell3)/dsqrt(mass((j-1)/3+1)*mass((i-1)/3+1))
!     ENDDO
!  ENDDO
!  GOTO 100
!200 CLOSE(50)
!
!  WRITE(*,*) 'End of ', TRIM(hessfile)
!
!!---- Build a dynamical matrix
!  DO k = -maxk,maxk
!     realk=dfloat(k)/dfloat(maxk)
!     DO i = 1, Ndof
!        DO j = 1, Ndof
!           df(i,j) = dcmplx(0.0d0,0.0d0)
!           DO m = -maxcell,maxcell
!              phase = dcmplx(0.0d0,-c_pi*dfloat(m)*realk)!check this part, sign doesn't change anything
!              df(i,j) = df(i,j) + f(i,j,kx*m,ky*m,kz*m) * cdexp(phase)
!           ENDDO !m
!        ENDDO !j
!     ENDDO !i
!     !-------- Diagonalize the matrix
!     nc=df
!     CALL complex_diag(Ndof,nc,freq)
!     !-------- Print the results
!     DO i = 1, Ndof
!        IF (freq(i).GE.0.0d0) THEN
!           freq(i) = dsqrt(freq(i))!*5140.4862d0
!        ELSE
!           freq(i) = -dsqrt(-freq(i))!*5140.4862d0
!        ENDIF
!
!        WRITE(10+i,*) realk,freq(i),dsqrt(f(1,1,0,0,0)/2.0)*2.0*dabs(dsin(0.5*c_pi*realk))
!        IF ((k.EQ.0).OR.(k.EQ.maxk)) THEN
!           x=0.0d0
!           y=0.0d0
!           z=0.0d0
!           DO j=1,Natom
!              x=x+cdabs(nc(3*j-2,i))**2
!              y=y+cdabs(nc(3*j-1,i))**2
!              z=z+cdabs(nc(3*j,i))**2
!           ENDDO !j
!           WRITE(*,'(i3,f10.2,a,3f8.5)') k,freq(i), " x,y,z =",x,y,z
!        ENDIF
!     ENDDO !i
!        k_freq(k,:)=freq(:)
!        k_coef(k,:,:)=nc(:,:)
!  ENDDO !k
!
!
!END SUBROUTINE dynamical
!
!SUBROUTINE xfc2qfc
!
!  USE constants, ONLY : c_pi
!  USE forceconstants
!  USE global
!  IMPLICIT NONE
!  INTEGER::k,k1,k2,k3,k4,m,m1,m2,m3,m4,i,i1,i2,i3,i4,c1,c2,c3,c4
!  COMPLEX*16 ::phase
!  ALLOCATE(x_grad(Ndof))
!  ALLOCATE(q_grad(Ndof,-maxk:maxk))
!
!  ALLOCATE(q_harm(Ndof,-maxk:maxk,Ndof,-maxk:maxk))
!  ALLOCATE(cck_coef(-maxk:maxk,Ndof,Ndof))
!          PRINT*,'allocated in xfc2qfc'
!  !  allocate(x_cube(Ndof,Ndof,-maxcell:maxcell,Ndof,-maxcell:maxcell))
!  !  allocate(q_cube(Ndof,-maxk:maxk,Ndof,-maxk:maxk,Ndof,-maxk:maxk))
!  !  allocate(x_quar(Ndof,Ndof,-maxcell:maxcell,Ndof,-maxcell:maxcell,Ndof,-maxcell:maxcell))
!  !  allocate(q_quar(Ndof,-maxk:maxk,Ndof,-maxk:maxk,Ndof,-maxk:maxk,Ndof,-maxk:maxk))
!  !  call read_fc()
!  q_grad=COMPLEX(0.0,0.0)
!  q_harm=COMPLEX(0.0,0.0)
!  cck_coef=CONJG(k_coef)
!    Ncell=maxcell*2+1
!    N_k=maxk*2+1
!  !  q_cube=complex(0.0,0.0)
!  !  q_quar=complex(0.0,0.0)
!  DO k=-maxk,maxk
!     IF (MOD(k,2).EQ. 0) THEN
!        DO m=1,Ndof
!           DO i=1,Ndof
!              q_grad(m,k)=q_grad(m,k)+ cck_coef(k,i,m)*x_grad(i)/dsqrt(mass(i))
!           ENDDO
!           q_grad(m,k)=q_grad(m,k)*dsqrt(dfloat(Ncell))
!        ENDDO
!     ENDIF
!  ENDDO
! PRINT*,'grad loop done'
!
!  DO k1=-maxk,maxk
!     DO k2=-maxk,maxk
!        IF (MOD(k1+k2,2*maxk).EQ. 0) THEN
!           DO m1=1,Ndof
!              DO m2=1,Ndof
!                 DO i1=1,Ndof
!                    DO i2=1,Ndof
!                       DO c1=-maxcell,maxcell
!                          phase=dcmplx(0.0,dfloat(c1*k2)*c_pi/dfloat(maxk))
!                          q_harm(m1,k1,m2,k2)=q_harm(m1,k1,m2,k2)+ &
!                               & cck_coef(k1,m1,i1)*cck_coef(k2,m2,i2)*x_harm(i1,i2,c1)*cdexp(phase)/dsqrt(mass(i1)*mass(i2))
!                       ENDDO
!                    ENDDO
!                 ENDDO
!                 PRINT*,m1, m2, q_harm(m1,k1,m2,k2),k_freq(k1,m1)**2.d0
!              ENDDO
!           ENDDO
!        ENDIF
!     ENDDO
!  ENDDO
!PRINT*,'harm loop done'
!END SUBROUTINE xfc2qfc

SUBROUTINE vscf
  USE global
  USE lapack
  USE yazar
  USE constants
  USE qVSCFintegrals
  IMPLICIT NONE

  REAL*8,ALLOCATABLE:: ham(:,:),energies(:)
  INTEGER:: mode,order,iter,maxiter
  REAL*8:: energy,oldenergy
  REAL*8,PARAMETER:: thresh=0.000001
  iter=1
  maxiter=20
  oldenergy=1000
  energy=0.0
  order=maxbasis(1)
  ALLOCATE(ham(order,order))
  ALLOCATE(energies(order))
  PRINT*,'const'
  WRITE(*,'(d50.25)')c_try
  WRITE(*,'(d50.25)')c_pi
  WRITE(*,'(d50.25)')c_gamma
  DO WHILE ((ABS(energy-oldenergy) .GT. thresh) .AND. iter .LT. maxiter )
     oldenergy=energy
     energy=0.0
     PRINT*,'**************************  iteration **************************',iter
     iter=iter+1
     DO mode=1,nDOF
        CALL meanFieldHamiltonian(mode,order,ham)

        PRINT*,'hamiltonian for mode ',mode
              CALL writematrix(order,ham)

           CALL diag(order,ham,energies)

          energy=energies(1)
      !  energy=dsqrt(utwo(mode))/2.d0+uzero(mode)
      !  hrmfreq(mode)=dsqrt(utwo(mode))
        PRINT*,"energy= ", energy!,energies(1)
          coef(mode,:)=ham(:,1)
        ! PRINT*,'eigvec'!ham(1,:)
        ! CALL writematrix(order,ham)
        !       DEALLOCATE(ham,energies)
     ENDDO
  ENDDO

END SUBROUTINE vscf

SUBROUTINE qvscf
  USE global
  USE lapack
  USE yazar
  USE constants
  USE qVSCFintegrals
  IMPLICIT NONE

  REAL*8,ALLOCATABLE:: ham(:,:),energies(:)
  INTEGER:: mode,order,iter,maxiter
  REAL*8:: energy,oldenergy
  REAL*8,PARAMETER:: thresh=0.000001
  iter=1
  maxiter=20
  oldenergy=1000
  energy=0.0
  order=maxbasis(1)
  ALLOCATE(ham(order,order))
  ALLOCATE(energies(order))
  PRINT*,'const'
  WRITE(*,'(d50.25)')c_try
  WRITE(*,'(d50.25)')c_pi
  WRITE(*,'(d50.25)')c_gamma
  DO WHILE ((ABS(energy-oldenergy) .GT. thresh) .AND. iter .LT. maxiter )
     oldenergy=energy
     energy=0.0
     PRINT*,'**************************  iteration **************************',iter
     iter=iter+1
     DO mode=1,nDOF
        CALL meanFieldHamiltonian(mode,order,ham)

        PRINT*,'hamiltonian for mode ',mode
        !      CALL writematrix(order,ham)
        PRINT*,utwo(mode),dsqrt(utwo(mode)),uzero(mode),kinetic(1,mode)

        energy=dsqrt(utwo(mode))/2.d0+ 4.d0*uzero(mode)
        hrmfreq(mode)=dsqrt(utwo(mode))
        PRINT*,"energy= ", energy!,energies(1)
        !  coef(mode,:)=ham(:,1)
        ! PRINT*,'eigvec'!ham(1,:)
        ! CALL writematrix(order,ham)
        !       DEALLOCATE(ham,energies)
     ENDDO
  ENDDO

END SUBROUTINE qvscf

SUBROUTINE meanFieldHamiltonian(mode,order,Hvscf)
  USE global
  USE forceconstants
  USE integral
  USE integral2
  IMPLICIT NONE
  INTEGER::mode,diff,i,j,m,m1,m2,order
  REAL*8::tmp,meanfieldpotential,meanpotkin
  REAL*8,DIMENSION(order,order):: Hvscf
  PRINT*, "forming VSCF Hamiltonian for mode" , mode
  meanpotkin=meanFieldPotential(mode)
  PRINT*,'meanpotkin',meanpotkin
  Hvscf=0.0
  DO i=1, maxbasis(mode)
     tmp=0.0
     tmp= tmp + (i - 0.5) * hrmfreq(mode) * 0.5;!kinetic term
     tmp= tmp + Hii(mode) * q_integral(i, mode, 2, 0);
     tmp= tmp + Qiiii(mode) * q_integral(i, mode, 4, 0);
     tmp= tmp + meanpotkin


     IF (nMR .GT. 1) THEN
        DO m=1,ndof
           IF (m .EQ. mode)CYCLE
           tmp= tmp + hij(m,mode) * q_integral(i, mode, 1, 0) * fullintegral(m, 1);
           tmp= tmp + qiijj(m,mode) * q_integral(i, mode, 2, 0) * fullintegral(m, 2);
           tmp= tmp + ciij(mode,m) * q_integral(i, mode, 2, 0) * fullintegral(m, 1);
           tmp= tmp + ciij(m,mode) * q_integral(i, mode, 1, 0) * fullintegral(m, 2);
           tmp= tmp + qiiij(mode,m) * q_integral(i, mode, 3, 0) * fullintegral(m, 1);
           tmp= tmp + qiiij(m,mode) * q_integral(i, mode, 1, 0) * fullintegral(m, 3);
        ENDDO
        IF (nMR .GT. 2) THEN

           DO m1=1,ndof-2
              IF (m1 .EQ. mode)CYCLE
              DO m2=m1+1,ndof-1
                 IF (m2 .EQ. mode) CYCLE
                 tmp= tmp + cijk(mode,m1,m2) * q_integral(i, mode, 1, 0) * fullintegral(m1, 1)  * fullintegral(m2, 1);
                 tmp= tmp + qiijk(mode,m1,m2) * q_integral(i, mode, 2, 0) * fullintegral(m1, 1) * fullintegral(m2, 1);
                 tmp= tmp + qiijk(m1,mode,m2) * q_integral(i, mode, 1, 0) * fullintegral(m1, 2)* fullintegral(m2, 1);
                 tmp= tmp + qiijk(m2,mode,m1) * q_integral(i, mode, 1, 0) * fullintegral(m1, 1) * fullintegral(m2, 2);
              ENDDO!!m2
           ENDDO!m1
        ENDIF!nMR>2
     ENDIF!nMR>1
     Hvscf(i,i)=tmp
     DO j=i+1,maxbasis(mode)!
        tmp=0.0
        diff = j - i;
        IF (diff > 4) EXIT !break;
        IF (diff .EQ. 2) tmp= tmp-SQRT(float((j-1) * (j - 2))) * hrmfreq(mode) * 0.25 !!kinetic term
        tmp= tmp + q_integral(j, mode, 1, diff) * gi(mode);
        tmp= tmp + q_integral(j, mode, 2, diff) * hii(mode);
        tmp= tmp + q_integral(j, mode, 3, diff) * ciii(mode);
        tmp= tmp + q_integral(j, mode, 4, diff) * qiiii(mode);
        ! print*,'test2',tmp
        IF (nMR .GT. 1)THEN
           DO m=1,ndof
              IF (m .EQ. mode)CYCLE
              tmp= tmp + hij(m,mode) * q_integral(j, mode, 1, diff) * fullintegral(m, 1);
              tmp= tmp + qiijj(m,mode) * q_integral(j, mode, 2, diff) * fullintegral(m, 2);
              tmp= tmp + ciij(mode,m) * q_integral(j, mode, 2, diff) * fullintegral(m, 1);
              tmp= tmp + ciij(m,mode) * q_integral(j, mode, 1, diff) * fullintegral(m, 2);
              tmp= tmp + qiiij(mode,m) * q_integral(j, mode, 3, diff) * fullintegral(m, 1);
              tmp= tmp + qiiij(m,mode) * q_integral(j, mode, 1, diff) * fullintegral(m, 3);
           ENDDO
           IF (nMR .GT. 2) THEN
              DO m1=1,ndof
                 IF (m1 .EQ. mode)CYCLE
                 DO m2=m1 + 1,ndof
                    IF (m2 .EQ. mode)CYCLE
                    tmp= tmp + cijk(mode,m1,m2) * q_integral(j, mode, 1, diff) * fullintegral(m1, 1)* fullintegral(m2, 1);
                    tmp= tmp + qiijk(mode,m1,m2) * q_integral(j, mode, 2, diff) * fullintegral(m1, 1)* fullintegral(m2, 1);
                    tmp= tmp + qiijk(m1,mode,m2) * q_integral(j, mode, 1, diff) * fullintegral(m1, 2)* fullintegral(m2, 1);
                    tmp= tmp + qiijk(m2,mode,m1) * q_integral(j, mode, 1, diff) * fullintegral(m1, 1)* fullintegral(m2, 2);
                 ENDDO!!m2
              ENDDO!!m1
           ENDIF!nMR>2
        ENDIF!nMR>1
        ! print*,'test3',tmp
        Hvscf(i,j)=tmp
     ENDDO!!end for j
  ENDDO!!end for i
END SUBROUTINE meanFieldHamiltonian

REAL*8 FUNCTION meanfieldpotential(mode)
  USE global
  USE forceconstants
  USE integral
  USE integral2
  USE constants
  IMPLICIT NONE
  REAL*8 :: tmp,meanFieldKinetic
  INTEGER::m,m1,m2,m3,mode
  tmp = 0.0;
  WRITE(*,'(d50.25)')c_gamma
  DO m=1,Ndof
     IF (m .EQ. mode)CYCLE
     tmp =tmp+ fullintegral(m, 1) * gi(m);
     print*,tmp
     tmp =tmp+ fullintegral(m, 2) * hii(m);
     print*,tmp
     tmp =tmp+ fullintegral(m, 3) * ciii(m);
     print*,tmp
     tmp =tmp+ fullintegral(m, 4) * qiiii(m);
     print*,tmp
     tmp =tmp+ meanFieldKinetic(m);
     print*,tmp
  ENDDO
  IF (nMR .GT. 1)THEN
     ! print*, "**************************nMR >1 "
     DO m1=1,nDOF
        IF (m1 .EQ. mode)CYCLE
        DO m2=m1+1,nDOF
           IF (m2 .EQ. mode .OR. m1 .EQ. m2)CYCLE
           tmp =tmp+ fullintegral(m1, 1) * fullintegral(m2, 1) * hij(m1,m2);
           print*,tmp
           tmp =tmp+ fullintegral(m1, 2) * fullintegral(m2, 2) * qiijj(m1,m2);
           print*,tmp
           tmp =tmp+ fullintegral(m1, 2) * fullintegral(m2, 1) * ciij(m1,m2);
           print*,tmp
           tmp =tmp+ fullintegral(m1, 1) * fullintegral(m2, 2) * ciij(m2,m1);
           print*,tmp
           tmp =tmp+ fullintegral(m1, 3) * fullintegral(m2, 1) * qiiij(m1,m2);
           print*,tmp
           tmp =tmp+ fullintegral(m1, 1) * fullintegral(m2, 3) * qiiij(m2,m1);
        print*,tmp
        ENDDO!m2
     ENDDO!m1
     IF (nMR .GT. 2)THEN
        DO m1=1,nDOF
           IF (m1 .EQ. mode) CYCLE
           DO m2=m1+1,nDOF
              IF (m2 .EQ. mode) CYCLE
              DO m3=m2+1,nDOF
                 IF (m3 .EQ. mode) CYCLE
                 tmp =tmp+ fullintegral(m1, 1) * fullintegral(m2, 1) * fullintegral(m3, 1) * cijk(m1,m2,m3);
                 !PRINT*,'tmp1 for mode ',m,' ',tmp !debug
                 tmp =tmp+ fullintegral(m1, 2) * fullintegral(m2, 1) * fullintegral(m3, 1) * qiijk(m1,m2,m3);
                 !PRINT*,'tmp2 for mode ',m,' ',tmp !debug
                 tmp =tmp+ fullintegral(m1, 1) * fullintegral(m2, 2) * fullintegral(m3, 1) * qiijk(m2,m1,m3);
                 !PRINT*,'tmp3 for mode ',m,' ',tmp !debug
                 tmp =tmp+ fullintegral(m1, 1) * fullintegral(m2, 1) * fullintegral(m3, 2) * qiijk(m3,m1,m2);
                 !PRINT*,'tmp4 for mode ',m,' ',tmp !debug
              ENDDO!m3
           ENDDO!m2
        ENDDO!m1
     ENDIF!nMR>2
  ENDIF!nMR>1
  meanfieldpotential=tmp
  print*, tmp,"meanfield pot for" ,mode
END FUNCTION meanfieldpotential

REAL*8 FUNCTION meanfieldkinetic(mode)
  USE global
  USE forceconstants
  USE integral
  USE integral2
  IMPLICIT NONE
  REAL*8 :: tmp
  INTEGER::i,j,m,m1,m2,m3,m4,diff,mode
  tmp = 0.0;
  DO i=1, maxbasis(mode)
     DO j=1, maxbasis(mode)
        diff = j - i;
        IF (diff .EQ. 0) THEN
           tmp =tmp+ coef(mode,i) * coef(mode,i) * (i - 0.5) * hrmfreq(mode) * 0.5;
        ELSE IF (ABS(diff) .EQ. 2) THEN
           tmp =tmp+ coef(mode,i) * coef(mode,j) * (-SQRT(float((MAX(i, j)-1) * (MAX(i, j) - 2)))) * hrmfreq(mode) * 0.25;!kinetic term
        ENDIF
     ENDDO!j
  ENDDO!i
  meanfieldkinetic=tmp
  print*, tmp,"meanfield kinetic for" ,mode
END FUNCTION meanfieldkinetic

SUBROUTINE init_QFF
  USE global, ONLY : nDOF,nMR
  USE forceconstants
  IMPLICIT NONE
  ALLOCATE(Gi(Ndof),Hii(ndof),Ciii(ndof),Qiiii(ndof))
  Gi=0.d0;Hii=0.d0;Ciii=0.d0;Qiiii=0.d0
  IF (nMR .GT. 1) THEN
     ALLOCATE(Hij(ndof,ndof),Ciij(ndof,ndof),Qiijj(ndof,ndof),Qiiij(ndof,ndof))
     Hij=0;Ciij=0;Qiijj=0;Qiiij=0
     IF (nMR .GT. 2) THEN
        ALLOCATE(Cijk(ndof,ndof,ndof),Qiijk(ndof,ndof,ndof))
        Cijk=0;Qiijk=0
     ENDIF
  ENDIF
END SUBROUTINE init_QFF
SUBROUTINE read_QFF
  USE global, ONLY : geo,nDOF,nMR
  USE forceconstants
  IMPLICIT NONE
  INTEGER :: i,j,k,n,tmp


  OPEN(unit=7,file='001.hs',status='OLD')
  READ(7,*,END=666)!# Energy / hartree
  READ(7,'(3x,D18.12)',END=666) E0
  READ(7,*,END=666)!# Geometry / Angs amu1/2
  j=MOD(Ndof,3)
  IF(j==0) THEN
     n=Ndof/3
     DO i=1,n
        READ(7,*,END=666) (geo(3*(i-1)+j),j=1,3)
     END DO
  ELSE
     n=(Ndof-j)/3
     DO i=1,n
        READ(7,*,END=666) (geo(3*(i-1)+j),j=1,3)
     END DO
     READ(7,*,END=666) (geo(j),j=3*n+1,Ndof)
  ENDIF

  !   >> >> 1 MR << <<
  READ(7,'(5x,a)') title1MR
  READ(7,*,END=666)!# Gradient / hartree Angs^-1 amu^-1/2
  DO n=1,Ndof
     READ(7,100)i, Gi(n)
  END DO
  READ(7,*,END=666)!# Hessian(i,i) / hartree Angs^-2 amu^-1
  DO n=1,Ndof
     READ(7,100,END=666) i,Hii(n)
  END DO
  READ(7,*,END=666)!# Cubic(i,i,i) / hartree Angs^-3 amu^-3/2
  DO n=1,Ndof
     READ(7,100,END=666)i, Ciii(n)
  END DO
  READ(7,*,END=666)!# Quartic(i,i,i,i) / hartree Angs^-4 amu^-2
  DO n=1,Ndof
     READ(7,100,END=666)i, Qiiii(n)
  END DO
100 FORMAT(i4,D20.10)

  !   >> >> 2 MR << <<
  IF(nMR .GT. 1) THEN
     tmp = ndof * (ndof - 1) / 2

     READ(7,'(5x,a)',END=222) title2MR

     READ(7,*,END=222)!# Hessian(i,j) / hartree Angs^-2 amu^-1
     DO n=1,tmp
        READ(7,200,END=666)i,j, Hij(i,j)
        Hij(j,i)=Hij(i,j)
     END DO

     READ(7,*,END=666)!# Quartic(i,i,j,j) / hartree Angs^-4 amu^-2
     DO n=1,tmp
        READ(7,200,END=666)i,j, Qiijj(i,j)
        Qiijj(j,i)=Qiijj(i,j)
     END DO

     READ(7,*,END=666)!# Cubic(i,i,j) / hartree Angs^-3 amu^-3/2
     DO n=1,2*tmp
        READ(7,200,END=666)i,j, Ciij(i,j)
     END DO

     READ(7,*,END=666)!# Quartic(i,i,i,j) / hartree Angs^-4 amu^-2
     DO n=1,2*tmp
        READ(7,200,END=666)i,j, Qiiij(i,j)
     END DO
200  FORMAT(2i4,D20.10)

     !   >> >> 3 MR << <<
     IF(nMR .GT. 2) THEN
        tmp = ndof * (ndof - 1) * (ndof - 2) / 6

300     FORMAT(3i4,D20.10)
        READ(7,'(5x,a)',END=222) title3MR
        READ(7,*,END=333)!# Cubic(i,j,k) / hartree Angs^-3 amu^-3/2
        DO n=1,tmp
           READ(7,300,END=666)i,j,k, Cijk(i,j,k)
           Cijk(i,k,j)=Cijk(i,j,k)
           Cijk(j,i,k)=Cijk(i,j,k)
           Cijk(j,k,i)=Cijk(i,j,k)
           Cijk(k,i,j)=Cijk(i,j,k)
           Cijk(k,j,i)=Cijk(i,j,k)
        END DO

        READ(7,*)!# Quartic(i,i,j,k) / hartree Angs^-4 amu^-2
        DO n=1,3*tmp
           READ(7,300,END=666)i,j,k, Qiijk(i,j,k)
           Qiijk(i,k,j)=Qiijk(i,j,k)
        END DO
        !  nMR=3
     ENDIF
  ENDIF
  RETURN

222 CONTINUE
  !  nMR=1
  RETURN
333 CONTINUE
  !  nMR=2
  RETURN
666 STOP 'error in QFF.pot'
END SUBROUTINE read_QFF

SUBROUTINE write_QFF
  USE global, ONLY : geo,nDOF
  USE forceconstants
  IMPLICIT NONE
  INTEGER :: i,j,k,n,tmp
  OPEN(unit=8,file='fc.poly')
  WRITE(8,1000)
  WRITE(8,'(3x,D18.12)') E0
  WRITE(8,1010)!# Geometry / Angs amu1/2
  j=MOD(Ndof,3)
  IF(j==0) THEN
     n=Ndof/3
     DO i=1,n
        WRITE(8,999) (geo(3*(i-1)+j),j=1,3)
     END DO
  ELSE
     n=(Ndof-j)/3
     DO i=1,n
        WRITE(8,999) (geo(3*(i-1)+j),j=1,3)
     END DO
     WRITE(8,999) (geo(j),j=3*n+1,Ndof)
  ENDIF
999 FORMAT(3d18.7)
1000 FORMAT('# Energy / hartree ')
1001 FORMAT('# Dipole / debye ')
1010 FORMAT('# Geometry / Angs amu1/2')

  !   >> >> 1 MR << <<

  WRITE(8,101) title1MR

  WRITE(8,110)!# Gradient / hartree Angs^-1 amu^-1/2
  DO n=1,Ndof
     WRITE(8,100)n, Gi(n)
  END DO

  WRITE(8,120)!# Hessian(i,i) / hartree Angs^-2 amu^-1
  DO n=1,Ndof
     WRITE(8,100)n, Hii(n)
  END DO

  WRITE(8,130)!# Cubic(i,i,i) / hartree Angs^-3 amu^-3/2
  DO n=1,Ndof
     WRITE(8,100)n, Ciii(n)
  END DO

  WRITE(8,140)!# Quartic(i,i,i,i) / hartree Angs^-4 amu^-2
  DO n=1,Ndof
     WRITE(8,100)n, Qiiii(n)
  END DO

100 FORMAT(i4,D20.10)
101 FORMAT('# 1MR ',a)
110 FORMAT('# Gradient / hartree Angs^-1 amu^-1/2 ')
111 FORMAT('# Gradient / debye Angs^-1 amu^-1/2 ')
120 FORMAT('# Hessian(i,i) / hartree Angs^-2 amu^-1 ')
121 FORMAT('# Hessian(i,i) / debye Angs^-2 amu^-1 ')
130 FORMAT('# Cubic(i,i,i) / hartree Angs^-3 amu^-3/2 ')
131 FORMAT('# Cubic(i,i,i) / debye Angs^-3 amu^-3/2 ')
140 FORMAT('# Quartic(i,i,i,i) / hartree Angs^-4 amu^-2 ')
141 FORMAT('# Quartic(i,i,i,i) / debye Angs^-4 amu^-2 ')

  !   >> >> 2 MR << <<

  tmp = ndof * (ndof - 1) / 2
  WRITE(8,201) title2MR

  WRITE(8,210)!# Hessian(i,j) / hartree Angs^-2 amu^-1
  DO i=2,ndof
     DO j=1,i-1
        WRITE(8,200)i,j, Hij(i,j)
     ENDDO
  END DO

  WRITE(8,220)!# Quartic(i,i,j,j) / hartree Angs^-4 amu^-2
  DO i=2,ndof
     DO j=1,i-1
        WRITE(8,200)i,j, Qiijj(i,j)
     END DO
  END DO

  WRITE(8,230)!# Cubic(i,i,j) / hartree Angs^-3 amu^-3/2
  DO i=2,ndof
     DO j=1,i-1
        WRITE(8,200)i,j, Ciij(i,j)
        WRITE(8,200)j,i, Ciij(j,i)
     END DO
  END DO

  WRITE(8,240)!# Quartic(i,i,i,j) / hartree Angs^-4 amu^-2
  DO i=2,ndof
     DO j=1,i-1
        WRITE(8,200)i,j, Qiiij(i,j)
        WRITE(8,200)j,i, Qiiij(j,i)
     END DO
  END DO

200 FORMAT(2i4,D20.10)
201 FORMAT('# 2MR ',a)
210 FORMAT('# Hessian(i,j) / hartree Angs^-2 amu^-1 ')
211 FORMAT('# Hessian(i,j) / debye Angs^-2 amu^-1 ')
220 FORMAT('# Quartic(i,i,j,j) / hartree Angs^-4 amu^-2 ')
221 FORMAT('# Quartic(i,i,j,j) / debye Angs^-4 amu^-2 ')
230 FORMAT('# Cubic(i,i,j) / hartree Angs^-3 amu^-3/2 ')
231 FORMAT('# Cubic(i,i,j) / debye Angs^-3 amu^-3/2 ')
240 FORMAT('# Quartic(i,i,i,j) / hartree Angs^-4 amu^-2 ')
241 FORMAT('# Quartic(i,i,i,j) / debye Angs^-4 amu^-2 ')
250 FORMAT(2i4,D20.10)

  !   >> >> 3 MR << <<
  tmp = ndof * (ndof - 1) * (ndof - 2) / 6
  WRITE(8,301) title3MR

  WRITE(8,310)!# Cubic(i,j,k) / hartree Angs^-3 amu^-3/2
  DO i=3,Ndof
     DO j=2,i-1
        DO k=1,j-1
           WRITE(8,300)i,j,k, Cijk(i,j,k)
        END DO
     END DO
  END DO

  WRITE(8,320)!# Quartic(i,i,j,k) / hartree Angs^-4 amu^-2
  DO i=3,Ndof
     DO j=2,i-1
        DO k=1,j-1
           WRITE(8,300)i,j,k, Qiijk(i,j,k)
           WRITE(8,300)j,k,i, Qiijk(j,k,i)
           WRITE(8,300)k,i,j, Qiijk(k,i,j)
        END DO
     END DO
  END DO

300 FORMAT(3i4,D20.10)
301 FORMAT('# 3MR ',a)
310 FORMAT('# Cubic(i,j,k) / hartree Angs^-3 amu^-3/2 ')
311 FORMAT('# Cubic(i,j,k) / debye Angs^-3 amu^-3/2 ')
320 FORMAT('# Quartic(i,i,j,k) / hartree Angs^-4 amu^-2 ')
321 FORMAT('# Quartic(i,i,j,k) / debye Angs^-4 amu^-2 ')

END SUBROUTINE write_QFF

SUBROUTINE make_quartic
  USE forceconstants
  USE global, ONLY: nMR
  IMPLICIT NONE
  Gi=0.0
  Ciii=0.0
  Qiiii=0.0
  IF(nMR .GT. 1) THEN
     Hij=0.0
     Ciij=0.0
     Qiiij=0.0
     IF(nMR .GT. 2) THEN
        Cijk=0.0
        Qiijk=0.0
     ENDIF
  ENDIF
END SUBROUTINE make_quartic

SUBROUTINE make_harmonic
  USE forceconstants
  USE global, ONLY: nMR
  IMPLICIT NONE
  Gi=0.0
  Ciii=0.0
  Qiiii=0.0
  IF(nMR .GT. 1) THEN
     Hij=0.0
     Ciij=0.0
     Qiiij=0.0
     Qiijj=0.0
     IF(nMR .GT. 2) THEN
        Cijk=0.0
        Qiijk=0.0
     ENDIF
  ENDIF
END SUBROUTINE make_harmonic

SUBROUTINE make_zero
  USE forceconstants
  USE global, ONLY: nMR
  IMPLICIT NONE
  Gi=0.0
  Hii=0.0
  Ciii=0.0
  Qiiii=0.0
  IF(nMR .GT. 1) THEN
     Hij=0.0
     Ciij=0.0
     Qiiij=0.0
     Qiijj=0.0
     IF(nMR .GT. 2) THEN
        Cijk=0.0
        Qiijk=0.0
     ENDIF
  ENDIF
END SUBROUTINE make_zero

SUBROUTINE make_units
  USE constants, ONLY: c_h2wn
  USE forceconstants
  USE global, ONLY: nMR
  IMPLICIT NONE
  Gi=Gi * c_h2wn
  Hii=Hii * 0.5 * c_h2wn
  Ciii=Ciii * c_h2wn / 6.0
  Qiiii=Qiiii * c_h2wn / 24.0
  IF(nMR .GT. 1) THEN
     Hij=Hij * c_h2wn
     Ciij=Ciij * c_h2wn * 0.5
     Qiiij=Qiiij * c_h2wn / 6.0
     Qiijj=Qiijj * c_h2wn * 0.25
     IF(nMR .GT. 2) THEN
        Cijk=Cijk * c_h2wn
        Qiijk=Qiijk * c_h2wn * 0.5
     ENDIF
  ENDIF
END SUBROUTINE make_units

SUBROUTINE xfc2qfc2

  USE constants, ONLY : c_pi
  USE global
  USE forceconstants
  IMPLICIT NONE
  INTEGER::k,k1,k2,k3,k4,m,m1,m2,m3,m4,i,i1,i2,i3,i4,c1,c2,c3,c4
  COMPLEX*16 ::phase
  ALLOCATE(x_grad(Ndof))
  ALLOCATE(q_grad(Ndof,-maxk:maxk))
  ALLOCATE(q_harm(Ndof,-maxk:maxk,Ndof,-maxk:maxk))
  ALLOCATE(cck_coef(-maxk:maxk,Ndof,Ndof))
  PRINT*,'allocated in xfc2qfc'
  !  allocate(x_cube(Ndof,Ndof,-maxcell:maxcell,Ndof,-maxcell:maxcell))
  !  allocate(q_cube(Ndof,-maxk:maxk,Ndof,-maxk:maxk,Ndof,-maxk:maxk))
  !  allocate(x_quar(Ndof,Ndof,-maxcell:maxcell,Ndof,-maxcell:maxcell,Ndof,-maxcell:maxcell))
  !  allocate(q_quar(Ndof,-maxk:maxk,Ndof,-maxk:maxk,Ndof,-maxk:maxk,Ndof,-maxk:maxk))
  !  call read_fc()
  q_grad=COMPLEX(0.0,0.0)
  q_harm=COMPLEX(0.0,0.0)
  cck_coef=CONJG(k_coef)
  Ncell=maxcell*2+1
  N_k=maxk*2+1
  !  q_cube=complex(0.0,0.0)
  !  q_quar=complex(0.0,0.0)
  DO k=-maxk,maxk
     IF (MOD(k,2).EQ. 0) THEN
        DO m=1,Ndof
           DO i=1,Ndof
              q_grad(m,k)=q_grad(m,k)+ cck_coef(k,i,m)*x_grad(i)/dsqrt(mass(i))
           ENDDO
           q_grad(m,k)=q_grad(m,k)*dsqrt(dfloat(Ncell))
        ENDDO
     ENDIF
  ENDDO
  PRINT*,'grad loop done'

  DO k1=-maxk,maxk
     DO k2=-maxk,maxk
        IF (MOD(k1+k2,2*maxk).EQ. 0) THEN
           DO m1=1,Ndof
              DO m2=1,Ndof
                 DO i1=1,Ndof
                    DO i2=1,Ndof
                       DO c1=-maxcell,maxcell
                          phase=dcmplx(0.0,dfloat(c1*k2)*c_pi/dfloat(maxk))
                          q_harm(m1,k1,m2,k2)=q_harm(m1,k1,m2,k2)+ &
                               & cck_coef(k1,m1,i1)*cck_coef(k2,m2,i2)*x_harm(i1,i2,c1)*cdexp(phase)/dsqrt(mass(i1)*mass(i2))
                       ENDDO
                    ENDDO
                 ENDDO
                 PRINT*,m1, m2, q_harm(m1,k1,m2,k2),k_freq(k1,m1)**2.d0
              ENDDO
           ENDDO
        ENDIF
     ENDDO
  ENDDO
  PRINT*,'harm loop done'
END SUBROUTINE xfc2qfc2

